I only checked changes in 2D mlebabeclap_adotx_centroid. It does not seem correct to me. The grad_eb_of_phi_on_centroids_extdir function uses anrmx and anrmy as follows:

Real dphidn = rhs(1)*nrmx + rhs(2)*nrmy;
return dphidn;

If the return value is scaled by a constant, it can be re-scaled later. But in your change, the ratio of nrmx to nrmy has changed. So I don't see how that can be correct.

In grad_eb_of_phi_on_centroids_extdir, rhs is the difference in phi in the x and y directions, but it has not been divided by a physical dx quantity yet. In other words, I am treating rhs and the output of grad_eb_of_phi_on_centroids_extdir as

rhs(1) = $\Delta_{x} \phi$ = dx $\frac{\partial\phi}{\partial x}$
rhs(2) = $\Delta_{y} \phi$ = dy $\frac{\partial\phi}{\partial y}$
grad_eb_of_phi_on_centroids_extdir = nrmx * dx * $\frac{\partial\phi}{\partial x}$ + nrmy * dy * $\frac{\partial\phi}{\partial y}$

This is why I scaled anrmx and anrmy the way I did. With anrmx = (apxm - apxp) dhx and anrmy = (apym - apyp) dhy, feb evaluates to

feb = beb * ( (apxm - apxp) dhx dx $\frac{\partial\phi}{\partial x}$ + (apym - apyp) dhy dy $\frac{\partial\phi}{\partial y}$ )

Apply definitions of dhx and dhy:

feb = B * beb * ( $\frac{\mathrm{apxm} - \mathrm{apxp}}{\mathrm{dx}}$ $\frac{\partial\phi}{\partial x}$ + $\frac{\mathrm{apym} - \mathrm{apyp}}{\mathrm{dy}}$ $\frac{\partial\phi}{\partial y}$)

This handles both 2D and 3D, with 2D changed to be more consistent with 3D.

So the proposed change in 2D is not just more consistency, right?

Yes, I think you are right, it's more than just consistency. The extra dxi factor is necessary to actually get spatial derivatives. This effect only matters when dx =/= dy. My apologizes for not adding a test to this PR. I can convert to draft if that's more appropriate.

No problems. I thought the 2D was just a programming style change and I was focusing on whether the 2D behavior has changed.